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A method Is presented for solving turbulent flow problems on three-dimensional unstructured grids. 
Spatial discretization Is accomplished by a cell-centered finite- volume formulation using an accurate lin- 
ear reconstruction scheme and upwind flux differencing. Time is advanced by an implicit backward- 
E uler time-stepping scheme. Flow turbulence effects are modeled by the Spalart-Alimaras one-equation 
model, which Is coupled with a wall function to reduce the number of cells In the sublayer region of the 
boundary layer. A systematic assessment of the method is presented to devise guidelines for more strate- 
gic application of the technology to complex problems. The assessment includes the accuracy in predic- 
tions of skin-friction coefficient, law-of-the-wall behavior, and surface pressure for a flat-plate turbulent 
boundary layer, and for the ONERA M6 wing under a high Reynolds number, transonic, separated flow 
condition. 


Introduction 

Significant advancements are being made toward solving 
complex viscous flows on three-dimensional configurations 
using unstructured-grid methodology [1-8], While solving 
such flows on highly-stretched tetrahedral cells is consider- 
ably more difficult than on hexahedral cells, the primary ad- 
vantage is derived from the greatly reduced grid generation 
times. Ref. [9] has demonstrated that ‘viscous’ grids can be 
easily generated on complex shapes by the Advancing Front 
/ Advancing Layers methodology (AFM/ALM). It is antici- 
pated that in the near future, viscous tetrahedral grids will 
be generated on complex geometries in a matter of days, as 
are in viscid tetrahedral grids today. 

The viscous, tetrahedral-based unstructured flow solution 
methodology is maturing along two tracks; node-centered 
and cell-centered schemes, each with their relative merits. 
Node-centered schemes exploit an efficient edge-based data 
structure, and have demonstrated multigrid and parallel 
computer implementations [2,5], but generally require large 
tetrahedral grids. Cell-centered schemes exploit geometric 
features of tetrahedra for constructing accurate spatial re- 
construction schemes, and provide comparable accuracy 
with fewer tetrahedra, but have not been extended to multi- 
grid or parallel architectures, and have exhibited some limi- 
tations in solution stability. 

There is a need for systematic assessments of the accu- 
racy and behavior of the various schemes. The present work 
focuses on an assessment of the upwind, tetrahedral cell- 
centered finite- volume scheme of Ref. [10]. This method is 
extended herein to include the Spalart-Alimaras one-equa- 
tion turbulence model, and the coupling of that model with a 
wall function to reduce the number of cells in the sublayer 
region of the boundary layer. It is anticipated that the wall 
function approach may be applicable to 3D separated flows 
since flow is not stagnant along separation lines. 
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The assessments will be derived from the flat-plate 
boundary layer problem, and the ONERA M6 wing at a 
high Reynolds number, transonic, separated flow condition. 
Key issues will be addressed related to applying a tetrahe- 
dral based, cell-centered Navier-Stokes method to turbu- 
lent-flow problems. The objectives of the study are to: 

1. assess the accuracy of computing turbulent-flow 
pressure distributions and skin friction coeffi- 
cients with tetrahedral cells, 

2. investigate the accuracy and utility of a wall func- 
tion formulation for computing 3-D high Rey- 
nolds number, transonic, separated flow with tet- 
rahedral cells, 

3. establish guidelines for generating unstructured, 
tetrahedral ‘viscous* grids for solving turbulent 
flow problems accurately and efficiently, 

4. demonstrate a mesh sequencing strategy for accel- 
erating solution convergence. 

Governing Equations 

The fluid motion is governed by the time-dependent 
Navier-Stokes equations for an ideal gas which express the 
conservation of mass, momentum, and energy for a com- 
pressible Newtonian fluid in the absence of external forces. 
The equations are given below in integral form for a 
bounded domain Q with the boundary 9Q 


fJJMJ F (Q) fidS = J J G (Q) AdS (1) 
n an do 


where Q = [p, pu, pv, pw, e 0 ] 
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Laminar viscosity ji, is computed by Sutherland's law. With 
the ideal gas assumption, the normalized values for pressure 
and temperature can be expressed as 

P = (Y~ 1) ( e o~ (° 2 + y2 + w 2 ) j 
and 

T = yp/p 

where y is the ratio of specific heats and is prescribed as 1.4 
for air. 

Numerical Procedure 

A finite-volume discretization is applied to Eq. 1 which 
results in a consistent approximation to the conservation 
laws where the time rate of change of the state vector Q 
within the domain Q is balanced by the net fluxes of F and 
G across the boundary surface dCl . The spatial domain is 
divided into a finite number of tetrahedral cells, with each 
element serving as a computational cell. Thus, the dis- 
cretized solution to Eq. 1 results in a set of volume-averaged 
state variables Q which are in balance with the area-aver- 
aged fluxes (inviscid and viscous) across the cell faces. 

Inviscid Fluxes 

Inviscid flux quantities are computed across each ceil 
face using the Roe [11] flux-difference splitting approach 
(FDS), or the Van Leer [12] flux- vector splitting technique 
(FVS). Spatial discretization is accomplished by a novel cell 
reconstruction process, which is based on an analytical for- 
mulation for computing solution gradients within tetrahe- 
dral cells. 

Cell reconstruction scheme 

The higher-order reconstruction scheme, derived in Ref. 
[10] and illustrated in Fig. 1, is based on a Taylor series ex- 
pansion of the cell-averaged solution to the cell face. A key 


X xy = (V + Hi) (u y + v x) . T xr = 0* + tt t ) K + wj 

V * (P + P t ) (v z + w y ) 



The equations are nondimensionalized with freestream 
reference values for density and a speed of sound . 
Here n x , fi y , and ft z are Cartesian components of the exte- 
rior surface unit normal ft on the boundary 3Q . The Carte- 
sian velocity components are u, v, and w in the x, y, and z 
directions, respectively. The term e 0 is the total energy per 
unit volume. The Prandtl number, Pr, is assigned a value of 
0.72, and the turbulent Prandtl number, Pr, , the value 0.9. 


• - Cell-centroid 

■ - Cell- vertex (determined by weighted average) 
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where = 0, 1st order 

= 1, higher order 

Fig. 1 Reconstruction stendl for tetrahedral cell-centered 
scheme. 
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(3) 


component of the scheme is the reconstruction of surround- 
ing cell-averaged data to a common vertex or node by a 
weighted averaging procedure. Reference [13] proposed a 
scheme based on an inverse-distance weighted averaging of 
the primitive variables from the cell centroid to the cell ver- 
tices. While this approach has proven to be both accurate 
and robust through wide application to inviscid problems, it 
is not fully second-order accurate in space. It has been 
shown in Ref. [14] to be approximately 1. 85-order accurate. 

As development efforts progressed toward solving the 
Navier-Stokes equations on highly stretched tetrahedral 
grids, it became evident that the accuracy of the inverse-dis- 
tance averaging scheme was not adequate. Thus, a fully sec- 
ond-order accurate averaging procedure was implemented 
which is based on work by Holmes and Connell [4] and 
Rausch, et. al. [15]. The procedure is derived by solving a 
constrained minimization problem to determine weight fac- 
tors which satisfy Laplacian relationships presented in Ref. 
[10]. The algorithm reconstructs to machine accuracy the 
exact values of a linear function at a node from surrounding 
cell-centered function values on an arbitrary tetrahedral 
grid. Furthermore, the simple universal formula shown in 
Fig. 1 for expanding the cell-centered data to the cell faces 
also reconstructs the exact value of a linear function to the 
cell face. Thus, the entire spatial reconstruction scheme is 
termed second-order accurate, which has been verified by 
Mitchell [14]. 

There is, however, an unresolved shortcoming to the 
Laplacian-weighted averaging scheme. Each weight factor 
is assumed to vary by some small perturbation from 1. In 
order to achieve an exact reconstruction on highly stretched 
cells, these perturbations can actually become on the order 
of one, thus resulting in some negative weight factors. 
While it can be demonstrated that the computed weight fac- 
tors produce an exact linear reconstruction, those with nega- 
tive values violate the principle of positivity, with a detri- 
mental impact on stability during convergence [16,17]. It is, 
thus, necessary to clip the weight factors between 0 and 2, 
thereby losing some of the exactness of the linear recon- 
struction, but ensuring a more stable scheme. 

Viscous Fluxes 

The viscous fluxes G(Q) are approximated at the cell-face 
centroids by linear reconstruction which provides a continu- 
ous representation of the solution variables across the cell 
faces. The stencil, presented by Mitchell [14], utilizes the 
averaged solution quantity at the three vertices of a cell 
face, q.„ q. 2 , and q^j, and the cell-centered values of the 
two cells sharing the face, q c j and q^ where qKp,u,v,w,p). 
The derivatives for u, v, w, and T in Eq. 1, e.g. for u„ u r 
and u v are derived from a Cramer’s rule solution to 
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Time Integration 

The viscous computations are advanced to steady state by 
the implicit time integration algorithm of Ref. [18]. The 
scheme uses the linearized, backward Euler time differenc- 
ing approach to update the solution at each time step for the 
set of equations 


lA]*{AQ} n = {R}° 


where 


r ji * _ V dR u 
[A] ~ At BQ 


The linear system of equations are solved at each time 
step with a subiterative procedure where the tetrahedral 
cells are grouped into “colors” (different from face-color- 
ing) such that no two cells share a common face. 

Thus, the solution is computed by solving for all the un- 
knowns in a particular color by a point-Jacobi subiteration 
step before proceeding to the next color. Since the solution 
of the unknowns for each group can depend on those from 
previously computed groups, a Gauss-Seidel-like effect is 
realized. The method has the advantage of being completely 
vectorizable. 

Because of the number of operations required to invert a 
matrix depends on the matrix bandwidth, the left-hand side 
of the system of linear equations is evaluated with first- 
order differencing to reduce both required storage and com- 
puter time. Convergence of the subiterations is further ac- 
celerated by using Van Leer’s Flux Vector Splitting (FVS) 
on the left-hand side. Thus in the present study, first-order 
differencing and FVS are applied to the left-hand side, and 
higher order differencing and FDS to the right-hand side. 
The viscous Jacobian terms are included in the left-hand 
side of the equation. 

It is necessary to store [A ] D , which is a 5X5 matrix for 
each cell, thus, storage requirements are 180 words/cell for 
the implicit code. The code requires 84 |is/cell/cycle on a 
CRAY Y-MP, or 37 ps/cell/cycle on a CRAY C-90, with 20 
subiterations and higher-order differencing. For compari- 
son, the block-structured code CFL3D [19] requires approx- 
imately 50 words/cell and 12|is/cell/cycle on a Cray C-90. 
While there may be some room for further improvement in 
resource requirements of the unstructured code, such codes 
are typically more computer-intensive because of their gen- 
eralized data structure. The success of this new technology 
will hinge on reducing the time and expense of generating 
viscous grids. 

Convergence acceleration 

Convergence to the steady state solution is accelerated by 
sacrificing the time accuracy of the scheme, and advancing 
the equations at each mesh point in time by the maximum 
permissible time step in that region. Even with such a local 
time stepping strategy, experience with solving 3-D viscous 
problems with the present cell -centered scheme has shown 
that maximum Courant, Friedrichs, Lewy (CFL) number is 
limited to approximately 25. This limitation is a conse- 
quence of violating the principle of positivity in weighting 
factors, as noted an earlier section and in Refs. [16, 17]. 

The inherent stability limitation can be improved by scal- 
ing the CFL number according to the deviation of cell as- 
pect ratio from the ideal value of an isotropic tetrahedron. 
This enables the dominate flow field to evolve quickly with 
the higher CFL numbers, while restricting the more temper- 
amental ‘viscous’ cells. A relation has been derived 
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where AR = [9 (V.) V (S t _J 3 j /ar,*.,. AR jfc „ = 8/ (373) , 
and fto is a scale factor. Here, V c is the cell volume, and 
is the area of the largest face of the cell. 

The computations presented in this paper were performed 
with a scale factor of f&= 6. Thus, for a prescribed setting of 
CFL=150, the actual CFL number will be linearly scaled 
between 25 for the thinnest tetrahedral cell to 150 for the 
most isotropic cell. The ultimate benefit of this procedure 
was a factor-of-two reduction in required solution cycles 
and, hence, computer time. 

Turbulence Model 


Spalart-Allmaras 

Closure of the Reynolds stress is provided by the one- 
equation Spalart-Allmaras (S-A) turbulence model [20]. 
This model is derived “using empiricism and arguments of 
dimensional analysis, Galilean invariance, and selective de- 
pendence on the molecular viscosity”. The model solves a 
partial differential equation (PDE) over the entire field for a 
working variable, v , from which the eddy viscosity, Pt, can 
be extracted. The PDE is solved separately from the flow 
equations using the same backward Euler time integration 
scheme, which results in a loosely coupled system. The pro- 
duction and destruction terms have been modified as recom- 
mended in Ref. [20] to ensure positive eddy viscosity 
throughout the computation. 

On ‘no-slip’ surfaces, the dependent variable v is set to 
zero. For tangent-flow surfaces, a zero gradient of the vari- 
able is applied. Far field boundary conditions are applied by 
extrapolating v from the interior for outflow boundaries, 
and taken from the free stream for the inflow. 

The S-A model requires that the distance of each cell to 
the nearest wall be provided for the near- wall damping 
terms for cells which are in proximity to ‘viscous’ surfaces. 
These distances are determined prior to code execution for 
cells in the “viscous” layers and contribute to only a small 
portion of the overall overhead. 

Wall Function 

The S-A model has been coupled with a wall function for- 
mulation to reduce the need for grid-resolving the flow in 
the sublayer portion of a turbulent boundary layer. With this 
approach, the inner region of the boundary layer is modeled 
by an analytical function which is matched with the numeri- 
cal solution in the outer region. This has the advantage of 1) 
significantly reducing memory requirement by eliminating 
a large portion of cells normally required to resolve the sub- 
layer, and 2) improving overall convergence by removing 
the thinner, more highly stretched cells which add stiffness 
to the solution process. A similar approach was successfully 
demonstrated in Ref. [21] where a two equation k-e turbu- 
lence model was coupled with a wall function in a modified 
version of the present code. 

The present implementation of a wall function exploits 
the inherent “structure” present in viscous unstructured 
grids produced by the Advancing Layers Method [9]. As ev- 


ident in Fig. 2, the tetrahedral vertices or nodes are aligned 
along rays emanating from the surface. 



Fig. 2 Inherent “structure” of thin-layer tetrahedral grids. 

The selected wall function is a law-of-the-wall expression 
[22] derived by Spalding in 1961. With a single function, it 
models the inner laminar sublayer, a transition region, and 
the intermediate logarithmic layer of the turbulent boundary 
layer: 


r + 2 * 3-i 

♦ + -kB icu* | + (KU ) (KU ) 

n = u+e e -1 -ku - — ' — - 

2 6 


(5) 


where the nondimensionalized terms are 


VM J p w 


u. 


Here p w ,p», are the fluid density and laminar viscosity on 
the surface, respectively, and IV,I the velocity magnitude at 
an adjacent point located a normal distance n, away; u* is 
the friction velocity; k=0.4 and, B=J.5. 

A face-centered, “slip” velocity boundary condition is de- 
termined by a two-step process. First as illustrated in Fig. 2, 
p w and p*, are assigned values from a boundary node, and 
IV, I is defined by the reconstructed velocity magnitude at 
the first connected node which is located n, above the sur- 
face. Eq.5 is then solved by Newton-Raphson iteration for 
w*, which is assumed to apply at the boundary nodes. 

Next, the computed friction velocities, w*, from three 
nodes comprising a boundary face (shaded surface in Fig. 2) 
are averaged to establish a face -centered value, and standard 
face-centered flow boundary condition quantities are pre- 
scribed for p w and p*. The parameter, n, , is the normal dis- 
tance to the centroid of the boundary tetrahedral cell. With 
these values, Eq. 5 is solved once again by Newton-Raph- 
son iteration for the velocity magnitude, IV, I. A slip-veloc- 
ity boundary condition is defined by assigning the new IV, I 
to the boundary face, and multiplying it with direction co- 
sines extracted from a standard inviscid-type flow-tangency 
velocity vector. 

A wall boundary condition for turbulent viscosity, which 
is required by the S-A PDE equation, is computed from a re- 
lation presented in Ref. [22] 
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where u + =IV,l/u, from the boundary face. 

The present implementation considers no adjustment to 
adiabatic wall density, which is important to high speed 
flows. This effect will be included in future work. 


Results 

Results are presented below for the flat-plate boundary 
layer problem, and the ONERA M6 wing at high Reynolds 
number, transonic, separated flow conditions. The normal 
grid spacing across the boundary layer is prescribed by the 
exponential function 


Aiij = An, (1 + <s (1 +6) J-, ) i l (7) 


such that nj = n^ + Anj.,. 

The parameter An, is the spacing of the first node above 
the surface, while a and b are parameters which control the 
growth. An initial estimate of the normal point distribution, 
n*, etc., is determined by experimenting with parameter 
variations on an assumed l/7th law velocity profile. 

Flat-Plate Boundary Layer 

The flat-plate boundary layer solution is used to assess 
the accuracy of the wall function in predicting flat-plate tur- 
bulent skin friction. The computations were made on quasi- 
2D unstructured grids for M«=o.5, Re^ixio 6 . 

Grid 1 was generated by constructing a 49X12 H-topol- 
ogy structured grid with a normal spacing defined by 
An,=0.001L, a= 0.3, and fc=0.07 in Eq. 7, which yields 
roughly 5 nodes across the boundary layer at x/L=0.5 and 
an approximate n + at the first node of 80. The resulting 
upper domain boundary (k=12) is located at 0.22L. The 2D 
grid was stacked spanwise in 0.02L increments to form 
three planes resulting in a 3D structured dual-channel grid 
(49X3X12) of H-H topology. Each hexahedral cell was sub- 
divided into 2 prismatic cells, which were further subdi- 
vided into 3 tetrahedra each to form the 3D unstructured 
grid with 15,552 cells. The “flat plate” was defined by a co- 
sine clustering between the “structured” indices 
15 £ / £ 49 along the k-\ boundary with in viscid flow pre- 
scribed on the h=l boundary ahead of the plate. Boundary 
conditions of constant entropy and constant total enthalpy 
were prescribed on the inflow plane, while an extrapolation 
condition was applied to the upper and exit domain bound- 
aries. A constant freestream pressure was also imposed on 
the exit plane. 

A second grid was generated in a similar manner as the 
fust to explore the lower limits of grid coarseness on solu- 
tion accuracy. Grid 2 was constructed from a 49X6 H-topol- 
ogy with the Eq. 7 parameters of An^O.OOlL, a=2.0, and 
b=0.07. This resulted in a 3D channel grid (49X3X6) with 
2,880 cells, and an upper domain boundary (k=6) also at 
0.22L. 

Fig. 3 portrays the effect of normal grid density on the 
law-of-the-wall behavior at x/L=0.5, Re^lxlO 6 , for the two 
grids. The plotted nodal solutions were reconstructed from 
the surrounding tetrahedral cells using the weighted averag- 


ing procedure discussed in an earlier section. Note that the 
fust nodal value is matched with the log layer at approxi- 
mately n + =80 for both grids. Grid 1 has 5 nodes across the 
boundary layer, while Grid 2 has 3 nodes. 



+ 
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Fig. 3 Effect of grid density on law-of-the-wall behavior for 
flat-plate boundary layer flow, x/L=0.5 (M«=0.5, Re^lXlO*). 

The true test of the methodology is reflected in the skin 
friction coefficient in Fig. 4. Grid 1 displays excellent 
agreement over 0.2<x/L<1.0 with the theoretical coefficient 
for fully turbulent flow, C f = 0.0583 (Re x )~ 1/5 , which is 
based on the l/7th power law assumption. Grid 2 does not 



x/L 


Fig. 4 Effect of grid density on skin-friction coefficient for flat- 
plate boundary layer flow (M.^.5, Re L =2X10 4 ). 


exhibit the same level of agreement, but is remarkably close 
considering its extreme grid coarseness across the boundary 
layer in Fig. 3. 

Based on experience with structured- grid computations, 
one would expect to need between 15 and 40 cells to ade- 
quately resolve turbulent boundary-layer flow, Thus, the re- 
sults of Figs. 3 and 4 require further analysis. As noted ear- 
lier, each hexahedral cell is subdivided into 2 prismatic 
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cells, which are each subdivided further into 3 tetrahedral 
cells. For a cell-centered scheme, each tetrahedron func- 
tions as one computational cell. Thus, a cell-averaged solu- 
tion resolves the solution at three vertical positions within 
each prismatic cell, (for a total of six locations within the 
hexahedral cell). In contrast, a cell-centered structured-grid 
code [19], or prismatic unstructured code [23], would re- 
solve the solution at only one vertical position within their 
respective cell layers. Hence in Fig. 3, there are actually 3 
tetrahedral centroids between each plotted solution point 
which contribute to those points through the reconstruction 
process mentioned earlier. One can conclude from this dis- 
cussion that it is more correct to consider the boundary layer 
as being resolved by 15 cells in Grid 1, and 9 cells in Grid 2, 
rather than by 5 and 3 nodes, respectively. 

One final note; the spurious behavior in the computed 
skin friction present near the plate leading edge 
( 0 £ x < 0.2) in Fig. 4 may be due to some numerical anom- 
alies of the weighted averaging scheme at the stagnation 
point where an inviscid surface suddenly changes to a vis- 
cous surface. The author plans to revisit this anomaly at a 
later date. The principal interest for the present study is in 
the fully developed turbulent flow over the remaining re- 
gion of the plate. 

ONERA M6 Wing 

Grid generation 

Tetrahedral viscous and inviscid grids were generated for 
the ONERA M6 wing using the VGRIDns code [9]. The 
VGRIDns code is based on the advancing- front method 
(AFM) for generating triangular surface mesh and tetrahe- 
dral volume cells. 

The distribution of surface and field grid points is con- 
trolled by a ‘structured’ background grid [24]. This trans- 
parent grid consists of Cartesian mesh overlaying the entire 
domain upon which the user prescribes ‘point’ and ‘line’ 
sources to impose the desired spacing distribution. Parame- 
ters are available to control cell size, and the direction and 
intensity of spatial variation. Cells can be stretched aniso- 
tropically in directions of small gradients in order to reduce 
the overall grid size. A smooth variation of spacing is 
achieved throughout the computational domain by solving 
an elliptic partial differential equation on the Cartesian 
mesh. The approach is analogous to modeling heat diffusion 
from discrete heat sources in a conducting medium. 

Thin-layered tetrahedra are generated in the ‘viscous’ re- 
gions by the advancing-layers method (ALM), which is 
based entirely on a modified AFM. The grid is marched 
away from the surface along smoothed vectors with a user 
prescribed distribution function, e.g. Eq. 7. As the cell sizes 
increase, the Cartesian background grid provides for a 
smooth transition to the remaining grid which is generated 
by the conventional AFM. 

The marching process of the ALM produces prismatic- 
like layers of grid which are subdivided into 3 tetrahedra 
within each “prism”, as illustrated in Fig. 2. As with the 
flat-plate boundary layer grids, each ‘prismatic* base cell is 
resolved spatially by three computational cells for a tetrahe- 
dral cell-centered scheme. 


Test matrix 

Several tetrahedral grids, eight thin-layered and one con- 
ventional, were generated for the ONERA M6 wing (see 
Tables 1 and 2). The spatial sources for the background grid 
were prescribed to produce a coarse (6483 triangles) and a 
fine (8956 triangles) surface grid distribution on the wing 
(see Fig. 5), and remained unchanged thereafter. Anisotro- 



Flg. 5a Upper surface triangulation of ONERA M6 wing, 
coarse viscous grid, 6483 triangles on wing surface. 



Fig. 5b Upper surface trlangulatlon of ONERA M6 wing, fine 
viscous grid, 8956 triangles on wing surface. 
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pic stretching of the surface grids was applied in the span- 
wise direction to reduce the total number of required cells 
while maintaining good chordwise resolution. A typical off- 
body distribution of volume grid for the coarse mesh is indi- 
rectly reflected in Fig. 6 by the centerplane grid of the WF2- 



Fig. 6 Center plane triangulation of ONERA M 6 wing, coarse 
grid, WF2-6(C). 


6(C) configuration. Note the smooth transition from the lay- 
ered ‘viscous’ grid to the conventional inviscid grid. As is 
evident in Fig. 6, the grid has characteristics of a structured 
O-mesh, since clustering of cells in the wake region has not 
been applied. Wake clustering is a topic for future research. 


Table 1 . - Designations for ONERA M6 viscous tetrahedral 
grids 


Initial Spacing 
(Arij/c R00T ) x 10 

Number of cells (nodes) across boundary 
layer 

=12 cells 
(4 nodes) 

=18 cells 
(6 nodes) 

=30 cells 
(10 nodes) 

0.6 

— 

FV-8* 

— 

1.350 

— 

WF1-6 

— 

2.025 

WF2-4 

WF2-6(C,F) 

WF2-10 

4.050 

— 

WF4-6 

— 

6.075 

— 

WF6-6 

— 


^includes two additional points in the inner layer 


Table 2 . - Parameters and characteristics of viscous grids. 


Grid 

a 

(aeeBq.7) 

b 

(see Eq. 7) 

Number 

surface 

triangles 

Number of 
cells, 
N m u 

(N&ell) Viscous 
Itviacid 

Inviscid 

— 

— 

6483 

258,768 

— 

FV-8 

0.5 

0.07 

6483 

414,038 

1.60 1 

WF1-6 

0.95 

0.07 

6483 

356,093 

1.38 

WF2-4 

2.2 

0.00 

6483 

324,356 

1.25 

WF2-6 (C) 

0.8 

0.07 

6483 

356,472 

1.38 

WF2-10 

0.2 

0.07 

6483 

463,968 

1.79 

WF4-6 

0.56 

0.07 

6483 

359,268 

1.39 

WF6-6 

0.432 

0.07 

6483 

362,311 

1.40 

WF2-6(F) 

0.8 

0.07 

8956 

578.556 

— 


The designations “FV” and “WF” in Tables 1 and 2 de- 
note “full viscous” with grid resolved sublayer, and “wall 
function” with non-grid resolved sublayer, respectively. 


The numerical nomenclature, e.g. 2-6, provides a nominal 
indicator of the (initial spacing)-(number of nodes across 
boundary layer) at the 0.5 mean aerodynamic chord for a 
Re^iUxio 6 . 

The full viscous grid, FV-8 was designed to have approx- 
imately the same number of nodes in the outer layer of the 
boundary layer as the WF2-6 grid, i.e. six nodes (18 tetrahe- 
dral layers), plus two additional nodes in the sublayer, for a 
total of 8 nodes (24 tetrahedral layers). 

A conventional inviscid grid was generated from the 
same wing surface grid, and with the same spatial source 
distributions as the viscous grids, thus, serving as a refer- 
ence for measuring the additional cells requirements for vis- 
cous grids as shown in Table 2. Note that the viscous grids 
require from 25-percent for the WF2-4 with 324,356 cells 
(57,490 nodes) to 79-percent for the WF2-10 with 463,968 
cells (80,927 nodes) more tetrahedral cells than the standard 
inviscid unstructured grid. It is obvious from this table that 
grid size can become rather large if more cells are needed 
across the boundary layer. This factor highlights the strong 
need for techniques, such as a wall function, to keep the 
‘viscous’ overhead down to manageable levels. 

A structured-grid computation was repeated from Ref. 
[19] for comparison with the unstructured results. The grid 
consisted of a 193x49x33 C-O mesh (294,912 hexahedral 
cells) with a minimum normal spacing over the wing of 
0.000015c R oor- This spacing matches that of the centroid of 
the surface tetrahedral cells in the FV-8 grid. Ref. [19] re- 
ports that this initial spacing resulted in an average n* of 4 
over the wing for M"=0.84, ot^3.06°, Re mjkC =11.7xl0 6 . 

Solution convergence 

All turbulent flow computations in this study were per- 
formed at the flow conditions of M«,=0.8447, a=5.06°, and 
Re m ^=l 1. 7x10 s , which represents a high Reynolds number, 
transonic, separated- flow condition. A typical solution con- 
vergence is shown in Fig. 7. 


R 

r- 

* 



Iteration 


Iteration 


Fig. 7 Solution convergence history for WF2-6(C), 
(M „=0.8447, a=5.06*, R«W=11.7X10*) 


The reason for the leveling off of the residual curve at 3- 
orders of magnitude reduction is not fully understood, but 
may be due to an unsteady nature of complex flow separa- 
tion in the wing-tip region. Note that the lift coefficient sets 
up quickly, but it is necessary to run the solution longer to 
allow for the separated region to evolve fully. 

Resource requirements for the computations are pre- 
sented in Table 3. All of the viscous cases were run with 
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CFL numbers starting at 20 and ramping up to 150 over 20 
cycles. The computations include the cell aspect ratio based 
variable CFL scaling strategy discussed earlier with a fsi=6 
applied. The time for the FV-8 case is based on mesh se- 
quencing which will be described in a later section. Unre- 
solved difficulties were encountered while attempting to 
start the full-viscous case from frees tream initial conditions. 


Table 3 . • Resource requirements for unstructured cases. 


Grid 

Number of 
cycles 

Cray C-90 
time, hours 

Memory, 

MW 

ovisdd 

300 

0.6 

46 

FV-8 

900/425 

5.9* 

75 

WF1-6 

1500 

5.78 

64 

WF2-4 

900 

3.05 

59 

WF2-6(Q 

2000 

7.34 

64 

WF2-10 

2500 

12.38 

84 

WF4-6 

1500 

5.75 

65 

WF6-6 

1500 

5.88 

66 

WF2-6TF) 

2500 

15.53 

104 


•obtained with mesh sequencing 


The structured-grid computation was performed with 
CFL3D using a multigrid strategy. The solution was con- 
verged in 500 cycles, and required 16MW and 0.66 Cray C- 
90 hours. 

Comparison of methodobgies 
Fig. 8 displays, the surface flow patterns for the FV-8 



Fig. 8 Surface “oil-flow” patterns for ONERA M6 wing, FV-8 
grid (M„=0.8447, a=5.0<S', Re m . c =I1.7X10 < ). 


case, which reveals a substantial shock-induced separation 
on the outboard portion of the wing. These patterns were de- 
termined from the reconstructed velocities at the first node 
above the wing surface. The general pattern shown in Fig. 8 
is representative of that from all the WF-series wall-func- 
tion solutions. 

Fig. 9 portrays a comparison of longitudinal C p distribu- 
tions for unstructured fUll-viscous and wall function solu- 


tions, a full-viscous structured solution [19] obtained with 
the CFL3D code, and a reference unstructured in viscid re- 
sult The comparisons are presented at the four span stations 



Fig. 9 C p distributions from unstructured and structured 
grids for ONERA M6 wing. (M m =0^447, o= 5.06 # , 
Re mac =l 1.7X10*). 


denoted on Fig. 8. The unstructured viscous results are in 
the best agreement with the experimental data of Ref. [25]. 
Furthermore, the wall-function solution, WF2-6(C), is in 
good agreement with the full-viscous solution, FV-8. 

The structured result in Fig. 9, which also employs the 
Spalart-Allmaras turbulence model, generally predicts the 
shock location too far forward and misses the aft-chord 
pressures. However, Ref. [19] demonstrated a strong depen- 
dence of the flow solution on the selection of turbulence 
model. Better agreement with data was shown in Ref. [19] 
using other turbulence models. As a note of caution, the 
good agreement of the unstructured results for q^O.90 may 
be fortuitous since the flow structure at the wing tip in Fig. 
8 is extremely complex and may exceed the capability of 
the one-equation S-A turbulence model. 

Surface grid sensitivity 

Fig. 10 shows the effect of surface grid density on the 
chordwise C p distributions at four chord stations. The WF2- 
6(C) & (F) have identical initial grid spacings. The sensitiv- 
ity to surface grid is small at the ti=0.65 and 0.80 stations 
where the separation is somewhat well behaved. Differ- 
ences are much larger in proximity to the more complex 
flow region for t]>0.90. The fine grid is in better agreement 
with the structured-grid result at the latter two stations, 
which once again highlights the strong sensitivities of the 
flow in that region. Although sensitivities to surface grid 
can be large in the complex tip flow region, the parametric 
study of normal grid spacing was performed on the coarser 
surface grid in order to reduce the overall computational ex- 
pense. 
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10 ° 10 1 10 2 10 3 10 4 10 5 


Fig. 10 Effect of surface grid density on C p distributions for 
ONERA M6 wing. (M_=0.8447, a=5.06®, Re =11.7X10*). 


Normal grid density 

Fig. 1 1 depicts the effect of normal grid density (a hori- 
zontal cut across Table 1) on the chordwise C p distributions 
at four chord stations. The three solutions, which represent 



x/c x/c 

Fig. 11 Effect of normal grid density on Cp distributions for 
ONERA M6 wing. (M„=0.8447, a=5.06®, R«w=l 1.7X10*). 


Fig. 12 Effect of normal grid density on law-of-th e-wall behav- 
ior at x/c=0.5, 2y/b=0.15 for ONERA M6 wing. (M„=0^447, 
a=5.06*, Rem^l 1.7X10*). 

Effect of initial grid spacing 

Fig. 13 shows the effect of initial grid spacing (a vertical 
cut through the WF-series of Table 1) on the chordwise C p 
distributions at four chord stations. Each grid is sized to 



x/c 


x/c 


12% 18% and 30-cells across the boundary layer, are gener- 
ally in good agreement with each other and the experimen- 
tal data. 

The law-of-the-wall behavior of the boundary layer for an 
attached-flow region of the wing (x/c=0.5, 2y/b=0.15), is 
plotted in Fig. 12. The fixed initial spacing yield n + =71 for 
all three cases, while there are 4, 6, and 10 nodes across the 
boundary layer corresponding to the WF2-4, -6, and -10, re- 
spectively. Recall that there are three tetrahedra between 
each nodal point contributing to the reconstruction of the 
solution to the nodes. 


Fig. 13 Effect of initial nodal spacing on C p distributions for 
ONERA M6 wing. (M«=0.8447, a=5.06 # , 11.7X10*). 


have approximately 6 nodes (18 tetrahedra) across the 
boundary layer at die midchord of the mean aerodynamic 
chord. For the test flow conditions, the initial grid spacings 
yield an n* of 48, 71, 146, and 218 for the first node of the 
WF1-6, 2-6, 4-6, and 6-6, respectively, at x/c=0.5, 2y / 
b=0.15. The sensitivity to initial spacing is negligible for all 
cases at the i]=0.65 and 0.80 stations, and for WF1-6 and 
WF2-6 at 11=0.90 and 0.95. As might be expected, the 
agreement with data deteriorates at the higher values of n + 
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for q>0.90. This result suggests that fairly large values of n + 
can be applied in conjunction with a wall function for more 
well behaved separated flows. However, more restricted 
values should be used in regions with complex 3D separated 
flow structures, such as spiral separations or primary saddle 
points, as depicted in Fig. 8. 

Mesh Sequencing 

The mesh sequencing strategy is often employed as a 
means of accelerating solution convergence. This strategy 
involves establishing the primary flow field relatively 
quickly around a configuration using a coarse mesh, then 
transferring that solution onto a finer mesh to complete the 
final grid-resolved solution. 

A demonstration of this procedure is provided for two of 
the more costly wall function solutions from Table 3, the 
WF2-6(F) and WF2-10, and for the full viscous case, the 
FV-8. Fig. 14 compares computer time requirements (in 


2 



-4 



Cray C-90 time, hrs 


Cray C-90 time, hrs 


Fig. 14 Effect of mesh sequencing on solution convergence for 
WF2-6(F) grid, (solid - single grid solution, dash - mesh 
sequenced solution), M oo =0JM47, 0=5.06°, Re mac = 11.7X10*. 


the order of 40- to 45-percent for the cases shown. (A sav- 
ings is not included for the FV-8 because of difficulties in 
obtaining a single-grid solution for that case.) An additional 
benefit is derived from the lower memory usage of the 
coarse-grid solution (59 megawords for the WF2-4), thus 
enabling primary flow to be setup more quickly while run- 
ning in smaller queues on heavily used computers. 

Concluding Remarks 

A systematic study has been initiated to assess the utiliza- 
tion of the cell-centered unstructured scheme for obtaining 
accurate solutions to the Navier-Stokes equations on three- 
dimensional configurations in an efficient manner. Closure 
to the flow equations is provided by a one-equation Spalart- 
Allmaras turbulence model, which is coupled with a wall 
function. 

Excellent accuracy in predicting the law-of-the-wall be- 
havior and surface skin friction coefficient with tetrahedral 
cells was demonstrated for the flat-plate boundary layer 
problem. The applicability of the tetrahedral-based wall 
function approach to 3D high Reynolds number, transonic, 
separated flow was validated on a parametric set of grids for 
the ONERA M6 wing. The validations were supported by 
comparisons with experimental data and a companion struc- 
tured-grid solution. The parametric study revealed that rea- 
sonable viscous solutions can be obtained with approxi- 
mately 25- to 80-percent more cells, hence memory, than a 
standard anisotropically stretched inviscid grid. Guidelines 
are established for prescribing an efficient distribution of 
normal grid spacing. A 40- to 45-percent solution conver- 
gence acceleration was demonstrated using a mesh sequenc- 
ing strategy. 

While the present study concludes with useful guidelines 
and better understanding of the base methodology, the next 
step of applying this knowledge to more complex geome- 
tries is important. Work is currently underway toward that 
end. 


Cray C-90 hrs) to obtain convergence of residual error and 
lift coefficient for the WF2-6(F) mesh. The solid curve ap- 
plies to the single-grid computation which took 15.5 Cray 
C-90 hours for 2500 cycles. The dashed line denotes the ap- 
plication of mesh sequencing, starting from the coarse grid 
WF2-4 solution at 900 cycles (see Table 3), interpolating 
that solution onto the WF2-6(F) grid, and continuing to run 
for another 700 cycles with CFL=150. The history plots in 
Fig. 14 do not reflect the additional computer time used for 
interpolating the solution from coarse to fine mesh. 

The full benefit of mesh sequencing is presented in Table 
4 for the three candidate cases, which includes the overhead 
of interpolating solutions. Note that the total savings is on 


Table 4 .- Resource requirements for mesh sequencing from 
coarse grid, WF2-4. 


Filler grid 

C-90 hrs for 
interpolating 
solution 

AddT 

Cycles 

beyond 

WF2-4 

Total 

solution time, 
Cray C-90 
hrs 

Percent 

savings 

WF2-6(F) 

2.1 

700 

9.5 

39 

WF2-10 

2.4 

300 

6.9 

44 

^V-8 

1.1 

425 

5.9 

— 
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